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INVERSE OF THE VANDERMONDE MATRIX WITH APPLICATIONS 

byL. Richard Turner 
Lewis Research Center 

SUMMARY 

The inverse of the Vandermonde matrix is given in the form of the product U“^L”^ 
of two triangular matrices by the display of generating formulas from which the elements 
of U”^ and L"^ may be directly computed. From these, it is directly possible to gen- 
erate special formulas for the approximation of linear transforms such as integrals, in- 
terpolants, and derivates for a variety of functions that behave as P(x)f(x), where P(x) 
is a polynomial and f(x) is generally not a poljmomial and may have localized singular- 
ities. 


INTRODUCTION 

The Vandermonde matrix, sometimes called an alternant matrix, comes from the 
approximation by a polynomial of degree n - 1 of a function f(x) with known values 
at n distinct values of the independent variable x. 

Many important fimctions of applied analysis cannot be well represented by polyno- 
mials, but these functions are accurately represented as products of polynomials by 
other fimctions that may not be analytic in the sense of function theory but can be effec- 
tively manipulated. 

To facilitate both the generation of working formulas for integration, interpolation, 
and differentiation and the calculation of other linear transformations, the Inverse of the 
Vandermonde matrix is presented in a simple matrix product form in which the elements 
are all computed by elementary explicit or recursive formulas. A few examples of ap- 
lications are given. 


ANALYSIS 


The Vandermonde matrix A arises as the matrix of coefficients required to 



evaluate the coefficients a. in any poljmomial approximation, as, for example, 

p n 1 

y = aQ + a^x + agX + . . . + a^^ ^x (1) 

to a fxmction y given at the n distinct points Xj^, X 2 , Xg, • • • , x^. The matrix A 
has the form 



If A“^ is known, the value of the coefficients is formally given by the expression 

[aQ, aj, & 2 , • • ^n-l] [^1’ ^2’ ^3’ ' ' *> yj 

where the brackets denote coliunn matrices and y. is equal to y(Xj). 

A simple form of the Inverse matrix A”^ is described in terms of the product 

where U“^ is an upper triangular matrix and L"^ is a lower triangular ma- 
trix. 

The Vandermonde matrix A has the determinant equal to TT (x. - x.) (ref. 1, p. 9) 

j>l ^ ^ 

and is nonsingular if all values of x. are distinct. It can, therefore, be factored into a 
lower triangular matrix L and an upper triangular matrix U where A is equal to LU. 
The factorization is imique if no row or column interchanges are made and if it is speci- 
fied that the diagonal elements of U are imity. 

The upper triangular factor U and the inverse L“ ^ of the lower triangular factor 
L are developed in reference 2, but the authors were content to depend on the evaluation 
of the numerical values for the coefficients a^ of equation (1) by the solutions of the 
equation 

u[aQ, a^, ag, • • •, = L ^[y^, yg, • • • , yj (4) 

in each case. 
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It is seen in reference 2 that, with a translation to matrix notation, the elements 1., 
-1 

of L are given by the relations 




£y=0 Ki 




‘.■TTr^ 


otherwise 


k= 




J 


(5) 


The explicit form of L“ ^ for a few rows and colmnns is 


L-l = 


Xj-Xg 


0 

1 


X£-Xi 


0 

0 


(Xj - X2>(Xj - Xg) (Xg - Xj)(Xg - Xg) (Xg - Xj)(Xg - Xg) 


( 6 ) 


It is asserted and proved herein that the elements Ujj of U“ ^ are given by the 


definition 


"N 


Uii=l 


Uii = 0 


'^ij=Vl,j-l-"i,j-l*j-l otherwise_ 


(7) 


where 


Uqj = 0 


The first few rows and columns of the asserted invp 0"' ^ are 
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1 


-XlXgXg 



X1X2 

-(Xj + X2) XjX2 + X2X3 + XgXj 

1 -(XJ+X2+X3) 

0 1 


( 8 ) 


ii* -i 

It is noted that the j column of U" does not depend on x^ but only on 
X, , x«, • • • , X. A proof that U”^ is as described by definition (7) is developed by 
showing that the product AU” is lower triangular and, therefore, equal to L. 

By direct computation, this is true for the Vandermonde matrix involving the two 
coordinates x^ and X 2 : 


1 

1 




(9) 


At this point, it can be assumed that, for the coordinates Xj^, 3 ^, Xg, • • • , x^^, the 
product of AU”^ is lower triangular, and the effect of adding the ordinate x^^^j and the 
corresponding rows and columns of A and U~ ^ can be considered. It is recalled that 
the n*"^ column of U'^ involves only the ordinates x^, X 2 , Xg, • • • , x^^ j, and that 

2 -1 

the inner products of the rows ( 1 , x^, x^ , • • • ) of A and column n of U are zero 
except for the diagonal element, which has the value (x^ - x^) (x^ - Xg) (x^ " Be- 
cause column n + 1 of U"^ defined by equation (7) is a linear combination of the ele- 
ments of column n, the inner products of the first n - 1 rows of A and the n + 1 col- 
umn of U"^ are all zero. It remains to show only that the element Pj^ of the prod- 
uct matrix vanishes. 

From the recvursive definition of columns of U~^, it is seen that the n + 1 column 
can be represented as the weighted sum of the two column matrices 

[0’"ln’ "2n» • • •’ "nn] - ""n^n’ "2n» "Sn’ ‘ ’ ’^nn’ The inner product of the 

column by the n*^ row of A, which is |l, x^, x^, • • • , xj^, produces two terms that 
are exactly equal but opposite in sign; hence, their sum is zero and p^^ = u^^ ^ 
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(of U”^) is zero. Therefore, the product AU~^ of the augmented matrix is lower trian- 
gular, which was to have been proved. 

PROPERTIES OF FACTOR MATRICES 

in general, the premultiplication by L”^ of a matrix of values of the function f(x), 
whose values are known at the ordinates [xj^, Xg, Xg, • • •] , produces the divided differ- 
ences used in Newton's interpolation formula (ref. 1, p. 3). 

It is also noted that the matrix does not involve the last ordinate. Because the 
ordinates at which the data Y are known may be arbitrarily identified, any sii^le 
ordinate may be left unspecified, and the elements of U"^ and all but the last row of 
L“ ^ can be computed independently of the location of the unspecified ordinate. 

It should also be mentioned that the ordinates x^^ may be complex numbers as long 
as their values are distinct. 

For any evenly spaced set of ordinates x^ where Xj+i = 1 + the matrix L"^ has 
the specific numerical form 



of which the elements have the values 


i 


ij 


(-1)^'^^ /i - 1\ 
(i-i)i Vi-V 


where 


i - 1 
3 - 1 


are the binomial coefficients. 
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For the set of ordinates Xj = i, the upper factor U"^ has the form 

1-1 2-6 
0 1 -3 11 

0 0 1-6 

0 0 0 1 


where, in this special case, 

where are the Stirling numbers of the first kind (ref. 3, p. 135). 

APPLICATIONS 

Because the Vandermonde matrix arises from the basic problem of passing a poly- 
nomial function through a given set of data, the results may be thought of as the approxi- 
mation of a function y(x) by a polynomial y(x). To the extent that y(x) is a reasonable 
approximation of y(x), it is possible to approximate linear transformations of y(x) by 
the corresponding transformation of y(x). 

Most of the classical formulas for numerical integration, interpolation, or differen- 
tiation can be generated directly. Exceptions are cases such as Gaussian and Chebysheff 
integrations, which require that the ordinates be especially selected by other considera- 
tions. 

Few examples of standard results will be displayed, but the techniques of their gen- 
eration should be clear from the special cases. Two cases are considered by which 
special formulas of high accuracy may be generated. 

Formulas for integration of products of functions and other related linear transforms 
may be developed as follows: Consider an integral or other linear transform of the prod- 
uct y(x)f(x), which is herein designated as T[y(x)f(x)]. The coefficients of the approxi- 
mation to y(x), namely, 

y(x) = + aj^x + agx^ + • • . 
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are given by the relation 


[aQ, a^, ag, • • •] = U'^L"^[y(x^), y(xg), y(x3), • • •] (12) 

Then, if it is possible to develop suitably computable ejq)ressions for T(x”f(x)), the 
transform T(y(x)f(x)) is given by 

(T(f(x)), T(xf(x)), . . .)u‘^L"^[yj, yg, yg, . . .] 


and, if y(x) is very nearly a polynomial, T(y(x)f(x)) is reasonably approximated by 
T(y(x)f(x)). 

Because the matrices U“^ and L"^ exist, an array of Lagrai^ian coefficients may 
be computed by the evaluation of 

(T(f(x)), T(xf(x)), T(x2f(x)), . . (13) 

independently of the actual values of y(Xj^). 

A second situation arises when y(x) is known, perhaps for analytical reasons, to be 
expressible by y(x) = p(x)f(x), where p(x) is a polynomial and f(x) is some known func- 
tion. For example, f(x) and x^f(x) may be Lesbegue integrable but cannot be accurately 
approximated by a polynomial. In this case, y(x) is approximated by 

y(x) = ^aQ + aj^x + agX^ + • • -)f(x) (14) 

and the coefficients a^^ are computable by the relation 



’y(xi) y(xg) 
f^’ f^’ 


(15) 


The desired Lagrangian coefficients to compute an approximate transform are generated 
by evaluating 

(T(f(x)), T(xf(x)), . . .) _1_, . . X (16) 

’ \f(Xl) f(X2) f(Xg) J 

where the braces denote a diagonal matrix. The resultii^ matrix of coefficients operates 
directly on the data 

[yp yg, y3» • • •] 
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NUMERICAL EXAMPLES 


(1) start with a classical example of generating the Adams -Moulton predictor cor- 
rector formulas, with f(x) = 1. For this example, the basis set x^ is [-3, -2, -1, o], 
and U~^ has the form 



and L"^ is given by definition (10). 
The predictor formula is 



The corrector formula is found from 



where h is th^ unit of step size. Both sets of coefficients are therefore generated by 
evaluating the expression 



1 

3 

1 

3 





It is convenient, at least for manual effort, to evaluate the expressions from left to r%ht. 
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in particular, 


h 



1 

2 

1 

2 




( 20 ) 



The leading term of the error series for these integrations is obtained by applying the 

41? If 

resulting operators to the data xf — and subtracting the correct integral for this 
5"*» ^ 24 

term, which is — . The data apart from the multiplicative constant are [81, 16, 1, 0], 
120 


251 5 ”** 

and the results are - — i h y for the predictor and 

720 

(2) As an example of a more general process, an 


19 5”’* 

— — h y for the corrector. 

720 

integral of the form 



y(0<i(l) 

V(i) 


should be considered, where y(x) is assumed to be regular and known and capable of 
being adequately approximated by a polynomial. This is an example of equation (14) with 

the data set x. = [o, 1, 2, • • • ] , with a scale h, and with f(x) = l/>^. The transform 


T(?"f(x)) = 



(^)(2n+l)/2 


Hence, the appropriate set of Lagrangian integration coefficients is given by 
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where U"^ has the form 



0 0 0 

1-12 
0 1-3 

0 0 1 



and L"^ has the form of definition (10). 

Because integrals for each ordinate x. are desired, the first matrix becomes 


CU“^ = h 



2 

2 

2 


3 

5 

7 

2V2 

iVi 

3 

5 

7 

2V3 

V 

2V3 

5 

^Vi . 

7 

/ 2 

2 


44 

/ 

3 

15 

105 

/ 2V2 

iVi 

3 

-ivi 

15 

i^Vi 

105 

\ 2V3 

2V3 

iVi 

5 

32 Vi 

35 


= C 


and the final result CU“^L"^ computed for three data points only is 



CU"^L"^ = JL 
15 



14 

leVi 

eVi 



( 25 ) 


Th e coefficients operate on the data y(0), y(h), and y(2h). 

(3) If the data for an integration are y(x), which is known to be of the form P(x)/^, 

and lim y(x)‘\^ = P(0) is known, a form may be derived to compute f y(?)d^, which 

can be obtained from example (2) by postmultiplying the general result of the form 
CU ^ by the diagonal matrix (l, 1, Vz, Vs, • • •). The resulting matrix operates 


on the data [y(0), y(h), y(2h), 


]. The matrix for three data points is 



18 

14 

-2V2 \ 

I2V2 

I6V2 

• 

I2V3 

eVi 

I2V6 / 

!, it is known that 

y(x) has the form 


(26) 


limit at zero of y(x)>^ is not known, similar formulas can be generated by considering 
the point set for x. of [h, 2h, 3h, • • • ]. In this case, U”^ has the form given in equa- 
tion (22), L ^ has the form of definition (10), the matrix C has the form of equation (23), 


and CU~ ^ has the form 


/ 


CU"^ = h 


\ 


2 

-A 

12 

-712 \ 



3 

5 

105 1 


2V2 

-2>^ 

3 

5 

-488^ 

105 

(27) 

2V3 

0 

5 

-393 Vi 1 

105 / 



For three points of data at h, 2h, and 3h, 
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CU'^L"^ = hi 



-56 

6 ' 

15 

15 

5 

^yT2 


iV^ 

15 

15 

5 

ilVi 

zlyTz 

iVi 


(28) 


This matrix is finally postmultiplied by the diagonal matrix |l, V 2 , Vij. The re- 
sult of this postmultiplication is 


h_ 

15 


' 68 

-56Vi 

lQ\lz 

52>6 

-68 

12yl6 

42 V 3 

-24Ve 

36 


(29) 


It is easy to generate weight coefficients for variotis processes of integration and 
differentiation or for the calculation of other transforms that are disadvantageous from 
the standpoint of error propagation. This may lead to computational instability if these 
weight coefficients are used within a larger problem, such as the nvunerical solution of 
a differential equation. 

While there are no simple rules for the analysis of the problem of stability, a crude 
measure of stability can be obtained by computing where w. de- 


notes the weighting coefficients of the numerical transform. 

As an example, the measure has been computed for each 

of the expository computations in this report. The expository computations are identi- 
fied in table I by equation and row. 

The mere magnitude of the stability measure is not an adequate guide to the effec- 
tiveness of a given quadrature formula but, in general, large values of the stability meas- 
ure are viewed imfavorably. 

It may be noted that the repeated use of the corrector formula of the Adams- 
Bashforth process (table I, row 2 of eq. (21)) leads to stability in the integration of sys- 
tems of differential equations if the step size is adequately small (ref. 4). Note that the 
stability factor for equation (21) is 1. 803. 
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TABLE I. - STABILITY FACTORS FOR 


SELECTED EXAMPLES 


■ 

Equation 

Row 

Stability factor, 

V“I <-"//( E-.)’' 

(21) 

1 

7.433 


2 

1.803 

(25) 

1 

1.322 


2 

1.156 


3 

1.039 

(26) 

1 

1.368 


2 

1.144 


3 

1.072 

(29) 

1 

9.443 


2 

5.177 


3 

3.473 


On the other hai\d, the alternate use of the open three-point Newton-Cotes formula 
with a stability factor of 1. 732 and Simpson's rule (closed three-point Newton-Cotes for- 
mula) with a stability factor of 1.225, as is done in Milne's schedule for the integration 
of systems of differential equations, is almost always imstable. 

In general, it will be observed that when strong inferences concerning the data such 
as the use of very many points in a simple quadrature formula or inferences concerning 
values of data outside the interval actually spanned by the data as in equation (29) or 
row 1 of equation (21) are made, the stability factors are large, and the resulting method 
may be a poor one. 


CONCLUSION 

The inverse A"^ of the Vandermonde matrix is given in the form of the product of 
two triangular matrices U”^L“^ by the display of generating formulas from which the 
elements of U”^ and may be directly computed. From these, it is directly pos- 

sible to generate special formulas for the approximation of linear transforms, such as 
integrals, interpolants, and derivates, for a variety of functions that behave as P(x)f(x), 
where P(x) is a polynomial and f(x) is generally not a polynomial and may be singular. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, May 9, 1966, 

129-04-06-02-22. 
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